Load in data

#load in global data set
globe_ice.df <- read.csv(file = "./global/data/global_sea_ice_1979-2022.csv", header = TRUE) %>% 
  as_tibble() %>% 
  rename("global_ice" = "Ice",
         "global_anomaly" = "Anomaly")

summary(globe_ice.df)
##       Date           Date2                Year          Month      
##  Min.   :197901   Length:518         Min.   :1979   Min.   : 1.00  
##  1st Qu.:198910   Class :character   1st Qu.:1989   1st Qu.: 3.00  
##  Median :200008   Mode  :character   Median :2000   Median : 6.00  
##  Mean   :200015                      Mean   :2000   Mean   : 6.48  
##  3rd Qu.:201105                      3rd Qu.:2011   3rd Qu.: 9.00  
##  Max.   :202202                      Max.   :2022   Max.   :12.00  
##    global_ice    global_anomaly 
##  Min.   :-9999   Min.   :-9999  
##  1st Qu.:   20   1st Qu.:   -1  
##  Median :   23   Median :    0  
##  Mean   :  -16   Mean   :  -39  
##  3rd Qu.:   25   3rd Qu.:    0  
##  Max.   :   28   Max.   :    2
#load in northern hemisphere data set
nohem_ice.df <- read.csv(file = "./global/data/nohem_sea_ice_1979-2022.csv", header = TRUE) %>% 
  as_tibble() %>% 
  rename("nohem_ice" = "Value",
         "nohem_anomaly" = "Anomaly")

summary(nohem_ice.df)
##       Date           Date2                Year          Month      
##  Min.   :197901   Length:519         Min.   :1979   Min.   : 1.00  
##  1st Qu.:198910   Class :character   1st Qu.:1989   1st Qu.: 3.00  
##  Median :200008   Mode  :character   Median :2000   Median : 6.00  
##  Mean   :200019                      Mean   :2000   Mean   : 6.47  
##  3rd Qu.:201106                      3rd Qu.:2011   3rd Qu.: 9.00  
##  Max.   :202203                      Max.   :2022   Max.   :12.00  
##    nohem_ice     nohem_anomaly  
##  Min.   :-9999   Min.   :-9999  
##  1st Qu.:    9   1st Qu.:   -1  
##  Median :   12   Median :    0  
##  Mean   :  -27   Mean   :  -39  
##  3rd Qu.:   14   3rd Qu.:    0  
##  Max.   :   16   Max.   :    1
#load in southern hemisphere data set
sohem_ice.df <- read.csv(file = "./global/data/sohem_sea_ice_1979-2022.csv", header = TRUE) %>% 
  as_tibble() %>% 
  rename("sohem_ice" = "Value",
         "sohem_anomaly" = "Anomaly")

summary(sohem_ice.df)
##       Date           Date2                Year          Month      
##  Min.   :197901   Length:519         Min.   :1979   Min.   : 1.00  
##  1st Qu.:198910   Class :character   1st Qu.:1989   1st Qu.: 3.00  
##  Median :200008   Mode  :character   Median :2000   Median : 6.00  
##  Mean   :200019                      Mean   :2000   Mean   : 6.47  
##  3rd Qu.:201106                      3rd Qu.:2011   3rd Qu.: 9.00  
##  Max.   :202203                      Max.   :2022   Max.   :12.00  
##    sohem_ice     sohem_anomaly  
##  Min.   :-9999   Min.   :-9999  
##  1st Qu.:    6   1st Qu.:    0  
##  Median :   12   Median :    0  
##  Mean   :  -27   Mean   :  -39  
##  3rd Qu.:   17   3rd Qu.:    0  
##  Max.   :   20   Max.   :    2
#load in hemispheric daily sea ice extent for duration calculations/plots
nohem_ice_daily.df <- read.csv(file = "./global/data/N_seaice_extent_daily_v3.0.csv", header = TRUE) %>% 
  as_tibble() %>% 
  rename("nohem_extent" = "Extent") %>% 
  mutate(Date = as.Date(Date)) %>% 
  filter(!Date < "1988-12-31" & Year != 2022)  #remove up to last day of 1987, years prior sampled every other day, need last day for lag calculation, also remove 2022 since incomplete data
  

summary(nohem_ice_daily.df)
##       Date                 Year          Month            Day      
##  Min.   :1988-12-31   Min.   :1988   Min.   : 1.00   Min.   : 1.0  
##  1st Qu.:1997-04-01   1st Qu.:1997   1st Qu.: 4.00   1st Qu.: 8.0  
##  Median :2005-07-01   Median :2005   Median : 7.00   Median :16.0  
##  Mean   :2005-07-01   Mean   :2005   Mean   : 6.52   Mean   :15.7  
##  3rd Qu.:2013-09-30   3rd Qu.:2013   3rd Qu.:10.00   3rd Qu.:23.0  
##  Max.   :2021-12-31   Max.   :2021   Max.   :12.00   Max.   :31.0  
##   nohem_extent  
##  Min.   : 3.34  
##  1st Qu.: 8.30  
##  Median :11.86  
##  Mean   :11.11  
##  3rd Qu.:14.09  
##  Max.   :16.25
sohem_ice_daily.df <- read.csv(file = "./global/data/S_seaice_extent_daily_v3.0.csv", header = TRUE) %>% 
  as_tibble() %>% 
  rename("sohem_extent" = "Extent") %>% 
  mutate(Date = as.Date(Date)) %>% 
  filter(!Date < "1988-12-31") %>%  #remove up to last day of 1987, years prior sampled every other day, need last day for lag calculation
  filter(!Date < "1988-12-31" & Year != 2022)  #remove up to last day of 1987, years prior sampled every other day, need last day for lag calculation, also remove 2022 since incomplete data

summary(sohem_ice_daily.df)
##       Date                 Year          Month            Day      
##  Min.   :1988-12-31   Min.   :1988   Min.   : 1.00   Min.   : 1.0  
##  1st Qu.:1997-04-01   1st Qu.:1997   1st Qu.: 4.00   1st Qu.: 8.0  
##  Median :2005-07-01   Median :2005   Median : 7.00   Median :16.0  
##  Mean   :2005-07-01   Mean   :2005   Mean   : 6.52   Mean   :15.7  
##  3rd Qu.:2013-09-30   3rd Qu.:2013   3rd Qu.:10.00   3rd Qu.:23.0  
##  Max.   :2021-12-31   Max.   :2021   Max.   :12.00   Max.   :31.0  
##   sohem_extent  
##  Min.   : 2.08  
##  1st Qu.: 6.09  
##  Median :12.59  
##  Mean   :11.68  
##  3rd Qu.:17.24  
##  Max.   :20.20
#load in sea level data
sea_level.df <- read.csv(file = "./global/data/epa-sea-level.csv", header = TRUE) %>% 
  as_tibble() %>% 
  rename("mean_sea_level" = "CSIRO.Adjusted.Sea.Level",
         "lowerCI_sea_level" = "Lower.Error.Bound",
         "upperCI_sea_level" = "Upper.Error.Bound")

summary(sea_level.df)
##      Date                Year      mean_sea_level lowerCI_sea_level
##  Length:35          Min.   :1979   Min.   :5.36   Min.   :5.10     
##  Class :character   1st Qu.:1988   1st Qu.:6.16   1st Qu.:5.91     
##  Mode  :character   Median :1996   Median :6.67   Median :6.39     
##                     Mean   :1996   Mean   :6.99   Mean   :6.72     
##                     3rd Qu.:2004   3rd Qu.:7.75   3rd Qu.:7.47     
##                     Max.   :2013   Max.   :9.33   Max.   :8.99     
##  upperCI_sea_level
##  Min.   :5.63     
##  1st Qu.:6.41     
##  Median :6.94     
##  Mean   :7.26     
##  3rd Qu.:8.03     
##  Max.   :9.66
#load in temp data
global_temp.df <- read.csv(file = "./global/data/monthly_globaltemp.csv", header = TRUE) %>% 
  as_tibble() %>% 
  arrange(-row_number()) %>%  #reverse row order so earliest years come first, matches sea ice df
  mutate(Date = as.Date(Date) %>% 
           as.yearmon()) #change to month year format
#join the three dataframes
total_ice.df <- inner_join(globe_ice.df, nohem_ice.df,
                           by = c("Date" = "Date",
                                  "Date2" = "Date2",
                                  "Year" = "Year",
                                  "Month" = "Month")) 

total_ice.df <- inner_join(total_ice.df, sohem_ice.df,
                           by = c("Date" = "Date",
                                  "Date2" = "Date2",
                                  "Year" = "Year",
                                  "Month" = "Month"))

#tidy up data
total_ice.df <- total_ice.df %>% 
  filter(global_ice != -9999) %>% #remove rows with -9999 (no value)
  dplyr::select(!Date) %>% 
  rename(Date = Date2) %>% 
  mutate(Date = str_replace_all(string = Date, 
                            pattern = "/", 
                            replacement = "-"))

total_ice.df <- total_ice.df %>% 
  mutate(Date = mdy(Date))
  

#check how summary looks now
summary(total_ice.df)
##       Date                 Year          Month         global_ice  
##  Min.   :1979-01-01   Min.   :1979   Min.   : 1.00   Min.   :16.3  
##  1st Qu.:1989-11-23   1st Qu.:1989   1st Qu.: 3.00   1st Qu.:20.3  
##  Median :2000-08-16   Median :2000   Median : 6.00   Median :23.5  
##  Mean   :2000-08-03   Mean   :2000   Mean   : 6.48   Mean   :23.0  
##  3rd Qu.:2011-05-08   3rd Qu.:2011   3rd Qu.: 9.00   3rd Qu.:25.3  
##  Max.   :2022-02-01   Max.   :2022   Max.   :12.00   Max.   :27.8  
##  global_anomaly     nohem_ice     nohem_anomaly      sohem_ice    
##  Min.   :-3.720   Min.   : 3.57   Min.   :-3.020   Min.   : 2.16  
##  1st Qu.:-0.725   1st Qu.: 8.56   1st Qu.:-0.843   1st Qu.: 5.82  
##  Median :-0.100   Median :12.10   Median :-0.185   Median :12.12  
##  Mean   :-0.279   Mean   :11.39   Mean   :-0.274   Mean   :11.57  
##  3rd Qu.: 0.390   3rd Qu.:14.30   3rd Qu.: 0.340   3rd Qu.:17.02  
##  Max.   : 1.650   Max.   :16.34   Max.   : 1.260   Max.   :19.76  
##  sohem_anomaly    
##  Min.   :-2.1300  
##  1st Qu.:-0.3200  
##  Median :-0.0350  
##  Mean   :-0.0059  
##  3rd Qu.: 0.3300  
##  Max.   : 1.8500
str(total_ice.df)
## tibble [516 × 9] (S3: tbl_df/tbl/data.frame)
##  $ Date          : Date[1:516], format: "1979-01-01" "1979-02-01" ...
##  $ Year          : int [1:516] 1979 1979 1979 1979 1979 1979 1979 1979 1979 1979 ...
##  $ Month         : int [1:516] 1 2 3 4 5 6 7 8 9 10 ...
##  $ global_ice    : num [1:516] 20.8 19.3 20.3 22.9 24.7 ...
##  $ global_anomaly: num [1:516] 1.39 0.95 0.88 1.4 1.24 1.6 1.4 0.82 0.35 0.13 ...
##  $ nohem_ice     : num [1:516] 15.4 16.2 16.3 15.4 13.9 ...
##  $ nohem_anomaly : num [1:516] 0.99 0.88 0.91 0.76 0.57 0.76 0.84 0.84 0.64 0.4 ...
##  $ sohem_ice     : num [1:516] 5.4 3.14 4 7.49 10.83 ...
##  $ sohem_anomaly : num [1:516] 0.4 0.07 -0.03 0.64 0.67 0.84 0.56 -0.02 -0.29 -0.28 ...

Modeling relationship between sea ice change and sea level rise

#prep sea ice df for comparison with sea level data
total_ice_march.df <- total_ice.df %>% 
  filter(Month == 3) %>% #only use ice data from March, matches collectiion of annual sea level data
  dplyr::select(Year, global_ice, global_anomaly) %>% 
  mutate(global_ice_change = global_ice - lag(global_ice)) %>% #create annual change variable
  filter(Year != 1979) #remove first year with missing change value

#remove Date column
sea_level.df <- sea_level.df %>% 
  dplyr::select(!Date)

#join dataframes
sea_levelandice.df <- inner_join(total_ice_march.df, sea_level.df,
                                 by = ("Year" = "Year"))

#check out summary
summary(sea_levelandice.df)
##       Year        global_ice   global_anomaly    global_ice_change
##  Min.   :1980   Min.   :17.6   Min.   :-1.8300   Min.   :-1.1500  
##  1st Qu.:1988   1st Qu.:19.1   1st Qu.:-0.3375   1st Qu.:-0.5875  
##  Median :1996   Median :19.4   Median :-0.0100   Median :-0.1950  
##  Mean   :1996   Mean   :19.4   Mean   :-0.0253   Mean   :-0.0085  
##  3rd Qu.:2005   3rd Qu.:19.9   3rd Qu.: 0.4525   3rd Qu.: 0.4750  
##  Max.   :2013   Max.   :20.5   Max.   : 1.0000   Max.   : 2.0900  
##  mean_sea_level lowerCI_sea_level upperCI_sea_level
##  Min.   :5.60   Min.   :5.34      Min.   :5.85     
##  1st Qu.:6.17   1st Qu.:5.92      1st Qu.:6.42     
##  Median :6.73   Median :6.46      Median :7.00     
##  Mean   :7.04   Mean   :6.76      Mean   :7.30     
##  3rd Qu.:7.75   3rd Qu.:7.48      3rd Qu.:8.04     
##  Max.   :9.33   Max.   :8.99      Max.   :9.66
#manual colors for legend
colors_icelevel <- c("Global Sea Level" = "#3399FF", "Global Sea Ice" = "#993399")

#plot sea ice and level together
sea_levelandice.df %>% 
  ggplot(aes(x = Year)) +
  geom_line(aes(y = global_ice, color = "Global Sea Ice")) +
  geom_line(aes(y = mean_sea_level*3, color = "Global Sea Level")) +
  geom_ribbon(aes(ymin = lowerCI_sea_level*3, ymax = upperCI_sea_level*3), fill = "#99CCFF", alpha = 0.25) +
  #add themes to clean up graph
  scale_y_continuous(
    # Features of the first axis
    name = "Global Sea Ice Extent (millions km^2)",
    # Add a second axis and specify its features
    sec.axis = sec_axis(trans = ~./3, name = "Mean Global Sea Level Change (in)")) +
  scale_color_manual("", #manual legend
                     values = colors_icelevel) + #set of colors matched with plot type
  theme(legend.position = "bottom",
        legend.key = element_rect(fill = "white"))

#construct simple linear regression models
mod6 <- lm(mean_sea_level ~ global_ice, data = sea_levelandice.df)
mod7 <- lm(mean_sea_level ~ global_ice_change + poly(global_ice_change, 2), data = sea_levelandice.df)

#model summary
summary(mod6)
## 
## Call:
## lm(formula = mean_sea_level ~ global_ice, data = sea_levelandice.df)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.484 -0.823 -0.221  0.589  2.437 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)   15.868      5.681    2.79   0.0087 **
## global_ice    -0.454      0.292   -1.56   0.1297   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.07 on 32 degrees of freedom
## Multiple R-squared:  0.0703, Adjusted R-squared:  0.0412 
## F-statistic: 2.42 on 1 and 32 DF,  p-value: 0.13
summary(mod7)
## 
## Call:
## lm(formula = mean_sea_level ~ global_ice_change + poly(global_ice_change, 
##     2), data = sea_levelandice.df)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -1.859 -0.656 -0.208  0.559  2.352 
## 
## Coefficients: (1 not defined because of singularities)
##                             Estimate Std. Error t value            Pr(>|t|)    
## (Intercept)                    7.037      0.170   41.31 <0.0000000000000002 ***
## global_ice_change              0.116      0.210    0.55              0.5839    
## poly(global_ice_change, 2)1       NA         NA      NA                  NA    
## poly(global_ice_change, 2)2    2.866      0.993    2.89              0.0071 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.993 on 31 degrees of freedom
## Multiple R-squared:  0.218,  Adjusted R-squared:  0.167 
## F-statistic: 4.32 on 2 and 31 DF,  p-value: 0.0222
#check residuals vs fitted
par(mfrow = c(2,2))
plot(mod6)

par(mfrow = c(2,2))
plot(mod7)

#plot relationships
sea_levelandice.df %>% 
  ggplot(aes(x = global_ice,  y = mean_sea_level)) +
  stat_smooth(formula = "y ~ poly(x,2)",
              method = "lm")

sea_levelandice.df %>% 
  ggplot(aes(x = global_ice_change,  y = mean_sea_level)) +
  stat_smooth(formula = "y ~ poly(x,2)",
              method = "lm")

Monthly sea ice and temperature

#convert date to month year to match temperature data
total_ice_yearmon.df <- total_ice.df %>% 
  mutate(Date = as.yearmon(Date))

#split  temp df into two dfsbased on source
global_temp_gis.df <- global_temp.df %>% 
  filter(Source %in% "GISTEMP")

global_temp_gcag.df <- global_temp.df %>% 
  filter(Source %in% "GCAG")

#join dataframes
ice_temp_gis.df <- inner_join(total_ice_yearmon.df, global_temp_gis.df,
                          by = ("Date" = "Date")) 

ice_temp_gcag.df <- inner_join(total_ice_yearmon.df, global_temp_gcag.df,
                          by = ("Date" = "Date"))

#manual colors for legend
colors_icetemp <- c("Global Temperature" = "#3399FF", "Global Sea Ice" = "#993399")

#plot sea ice and level together
ice_temp_gis.df %>% 
  ggplot(aes(x = Date)) +
  geom_line(aes(y = global_ice, color = "Global Sea Ice")) +
  geom_line(aes(y = Mean*25, color = "Global Temperature")) +
  #add themes to clean up graph
  scale_y_continuous(
    # Features of the first axis
    name = "Global Sea Ice Extent (millions km^2)",
    # Add a second axis and specify its features
    sec.axis = sec_axis(trans = ~./25, name = "Mean Global Temp Change (degrees Celsius)")) +
  scale_color_manual("", #manual legend
                     values = colors_icetemp) + #set of colors matched with plot type
  theme(legend.position = "bottom",
        legend.key = element_rect(fill = "white"))

 #plot sea ice and level together
ice_temp_gis.df %>% 
  ggplot(aes(x = Date)) +
  geom_line(aes(y = global_anomaly, color = "Global Sea Ice")) +
  geom_line(aes(y = Mean*3, color = "Global Temperature")) +
  #add themes to clean up graph
  scale_y_continuous(
    # Features of the first axis
    name = "Global Sea Ice Anomalies",
    # Add a second axis and specify its features
    sec.axis = sec_axis(trans = ~./3, name = "Mean Global Temp Change (degrees Celsius)")) +
  scale_color_manual("", #manual legend
                     values = colors_icetemp) + #set of colors matched with plot type
  theme(legend.position = "bottom",
        legend.key = element_rect(fill = "white"))                   

#construct simple linear regression models
mod8 <- lm(global_ice ~ Mean, data = ice_temp_gis.df)
mod9 <- lm(global_anomaly ~ Mean, data = ice_temp_gis.df)

#model summary
summary(mod8)
## 
## Call:
## lm(formula = global_ice ~ Mean, data = ice_temp_gis.df)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -6.227 -2.568  0.848  2.146  4.510 
## 
## Coefficients:
##             Estimate Std. Error t value             Pr(>|t|)    
## (Intercept)   24.383      0.287   85.01 < 0.0000000000000002 ***
## Mean          -2.571      0.546   -4.71            0.0000033 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 2.75 on 452 degrees of freedom
## Multiple R-squared:  0.0468, Adjusted R-squared:  0.0447 
## F-statistic: 22.2 on 1 and 452 DF,  p-value: 0.00000327
summary(mod9)
## 
## Call:
## lm(formula = global_anomaly ~ Mean, data = ice_temp_gis.df)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.8921 -0.3729  0.0605  0.4468  1.6615 
## 
## Coefficients:
##             Estimate Std. Error t value            Pr(>|t|)    
## (Intercept)   0.6942     0.0694    10.0 <0.0000000000000002 ***
## Mean         -1.6366     0.1319   -12.4 <0.0000000000000002 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.666 on 452 degrees of freedom
## Multiple R-squared:  0.254,  Adjusted R-squared:  0.252 
## F-statistic:  154 on 1 and 452 DF,  p-value: <0.0000000000000002
#check residuals vs fitted
par(mfrow = c(2,2))
plot(mod8)

par(mfrow = c(2,2))
plot(mod9)

#plot relationships
ice_temp_gis.df %>% 
  ggplot(aes(x = Mean,  y = global_ice)) +
  stat_smooth(formula = "y ~ x",
              method = "lm")

ice_temp_gis.df %>% 
  ggplot(aes(x = Mean,  y = global_anomaly)) +
  stat_smooth(formula = "y ~ x", 
              method = "lm")

Spatial plots of sea ice

Read in a few csvs to get started

#load shapefiles
#northern hemisphere
Nseaice.mar1980 <- st_read("./global/data/Monthly Sea Ice Extent/all_spatial_files/extent_N_198003_polygon_v3.0.shp")
## Reading layer `extent_N_198003_polygon_v3.0' from data source 
##   `/Users/jackrabe/Documents/PUBH 7462 Advanced Programming in R/Final Project/ice.trends.github.io/global/data/Monthly Sea Ice Extent/all_spatial_files/extent_N_198003_polygon_v3.0.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 122 features and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -3375000 ymin: -4900000 xmax: 3600000 ymax: 5850000
## Projected CRS: NSIDC Sea Ice Polar Stereographic North
Nseaice.mar2020 <- st_read("./global/data/Monthly Sea Ice Extent/all_spatial_files/extent_N_202003_polygon_v3.0.shp")
## Reading layer `extent_N_202003_polygon_v3.0' from data source 
##   `/Users/jackrabe/Documents/PUBH 7462 Advanced Programming in R/Final Project/ice.trends.github.io/global/data/Monthly Sea Ice Extent/all_spatial_files/extent_N_202003_polygon_v3.0.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 162 features and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -3400000 ymin: -4875000 xmax: 3600000 ymax: 5850000
## Projected CRS: NSIDC Sea Ice Polar Stereographic North
Nseaice.sep1980 <- st_read("./global/data/Monthly Sea Ice Extent/all_spatial_files/extent_N_198009_polygon_v3.0.shp")
## Reading layer `extent_N_198009_polygon_v3.0' from data source 
##   `/Users/jackrabe/Documents/PUBH 7462 Advanced Programming in R/Final Project/ice.trends.github.io/global/data/Monthly Sea Ice Extent/all_spatial_files/extent_N_198009_polygon_v3.0.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 118 features and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -2550000 ymin: -3250000 xmax: 2350000 ymax: 2100000
## Projected CRS: NSIDC Sea Ice Polar Stereographic North
Nseaice.sep2020 <- st_read("./global/data/Monthly Sea Ice Extent/all_spatial_files/extent_N_202009_polygon_v3.0.shp")
## Reading layer `extent_N_202009_polygon_v3.0' from data source 
##   `/Users/jackrabe/Documents/PUBH 7462 Advanced Programming in R/Final Project/ice.trends.github.io/global/data/Monthly Sea Ice Extent/all_spatial_files/extent_N_202009_polygon_v3.0.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 83 features and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -2525000 ymin: -2800000 xmax: 2350000 ymax: 4950000
## Projected CRS: NSIDC Sea Ice Polar Stereographic North
#southern hemisphere
Sseaice.mar1980 <- st_read("./global/data/Monthly Sea Ice Extent/all_spatial_files/extent_S_198003_polygon_v3.0.shp")
## Reading layer `extent_S_198003_polygon_v3.0' from data source 
##   `/Users/jackrabe/Documents/PUBH 7462 Advanced Programming in R/Final Project/ice.trends.github.io/global/data/Monthly Sea Ice Extent/all_spatial_files/extent_S_198003_polygon_v3.0.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 31 features and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -2525000 ymin: -2475000 xmax: 2625000 ymax: 2325000
## Projected CRS: NSIDC Sea Ice Polar Stereographic South
Sseaice.mar2020 <- st_read("./global/data/Monthly Sea Ice Extent/all_spatial_files/extent_S_202003_polygon_v3.0.shp")
## Reading layer `extent_S_202003_polygon_v3.0' from data source 
##   `/Users/jackrabe/Documents/PUBH 7462 Advanced Programming in R/Final Project/ice.trends.github.io/global/data/Monthly Sea Ice Extent/all_spatial_files/extent_S_202003_polygon_v3.0.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 12 features and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -2500000 ymin: -2600000 xmax: 2875000 ymax: 2300000
## Projected CRS: NSIDC Sea Ice Polar Stereographic South
Sseaice.sep1980 <- st_read("./global/data/Monthly Sea Ice Extent/all_spatial_files/extent_S_198009_polygon_v3.0.shp")
## Reading layer `extent_S_198009_polygon_v3.0' from data source 
##   `/Users/jackrabe/Documents/PUBH 7462 Advanced Programming in R/Final Project/ice.trends.github.io/global/data/Monthly Sea Ice Extent/all_spatial_files/extent_S_198009_polygon_v3.0.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 5 features and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -3050000 ymin: -3025000 xmax: 3825000 ymax: 4000000
## Projected CRS: NSIDC Sea Ice Polar Stereographic South
Sseaice.sep2020 <- st_read("./global/data/Monthly Sea Ice Extent/all_spatial_files/extent_S_202009_polygon_v3.0.shp")
## Reading layer `extent_S_202009_polygon_v3.0' from data source 
##   `/Users/jackrabe/Documents/PUBH 7462 Advanced Programming in R/Final Project/ice.trends.github.io/global/data/Monthly Sea Ice Extent/all_spatial_files/extent_S_202009_polygon_v3.0.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 3 features and 1 field
## Geometry type: POLYGON
## Dimension:     XY
## Bounding box:  xmin: -2825000 ymin: -3050000 xmax: 3350000 ymax: 4000000
## Projected CRS: NSIDC Sea Ice Polar Stereographic South
#make basic plot of layers
#northern hemisphere
#difference in ice between mar 1980 and mar 2020
 ggplot() +
  geom_sf(data = Nseaice.mar1980, fill = "#3399FF", show.legend = TRUE) +
  geom_sf(data = Nseaice.mar2020, fill = "#9900FF", show.legend = TRUE) +
  scale_fill_manual("Sea Ice Extent", values = c("Mar 1980" = "#3399FF", "Mar 2020" = "#9900FF"))

#difference in ice between sept 1980 and sept 2020
ggplot() +
  geom_sf(data = Nseaice.sep1980, fill = "#FF6600", show.legend = TRUE) +
  geom_sf(data = Nseaice.sep2020, fill = "#33CC33", show.legend = TRUE) +
  scale_fill_manual("Sea Ice Extent", values = c("Sept 1980" = "#FF6600", "Sept 2020" = "#33CC33"))

#difference in ice between mar 2020 and sept 2020
ggplot() +
  geom_sf(data = Nseaice.mar2020, fill = "#996633", show.legend = TRUE) +
  geom_sf(data = Nseaice.sep2020, fill = "#FFFF00", show.legend = TRUE) +
  scale_fill_manual("Sea Ice Extent", values = c("Mar 2020" = "#996633", "Sept 2020" = "#FFFF00"))

#outhern hemisphere
#difference in ice between mar 1980 and mar 2020
ggplot() +
  geom_sf(data = Sseaice.mar1980, fill = "#3399FF", show.legend = TRUE) +
  geom_sf(data = Sseaice.mar2020, fill = "#9900FF", show.legend = TRUE) +
  scale_fill_manual("Sea Ice Extent", values = c("Mar 1980" = "#3399FF", "Mar 2020" = "#9900FF"))

#difference in ice between sept 1980 and sept 2020
ggplot() +
  geom_sf(data = Sseaice.sep1980, fill = "#FF6600", show.legend = TRUE) +
  geom_sf(data = Sseaice.sep2020, fill = "#33CC33", show.legend = TRUE) +
  scale_fill_manual("Sea Ice Extent", values = c("Sept 1980" = "#FF6600", "Sept 2020" = "#33CC33"))

#difference in ice between mar 2020 and sept 2020
ggplot() +
  geom_sf(data = Sseaice.sep2020, fill = "#FFFF00", show.legend = TRUE) +
  geom_sf(data = Sseaice.mar2020, fill = "#996633", show.legend = TRUE) +
  scale_fill_manual("Sea Ice Extent", values = c("Mar 2020" = "#996633", "Sept 2020" = "#FFFF00"))

#load shapefile
seaice.jan1980 <- readOGR("./global/data/Monthly Sea Ice Extent/all_spatial_files/extent_N_198001_polygon_v3.0.shp")
## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/jackrabe/Documents/PUBH 7462 Advanced Programming in R/Final Project/ice.trends.github.io/global/data/Monthly Sea Ice Extent/all_spatial_files/extent_N_198001_polygon_v3.0.shp", layer: "extent_N_198001_polygon_v3.0"
## with 143 features
## It has 1 fields
## Integer64 fields read as strings:  FID
#check class and crs
class(seaice.jan1980)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
crs(seaice.jan1980)
## Coordinate Reference System:
## Deprecated Proj.4 representation:
##  +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +x_0=0 +y_0=0 +a=6378273
## +rf=298.279411123064 +units=m +no_defs +type=crs 
## WKT2 2019 representation:
## PROJCRS["NSIDC Sea Ice Polar Stereographic North",
##     BASEGEOGCRS["Unspecified datum based upon the Hughes 1980 ellipsoid",
##         DATUM["Not specified (based on Hughes 1980 ellipsoid)",
##             ELLIPSOID["Hughes 1980",6378273,298.279411123064,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4054]],
##     CONVERSION["US NSIDC Sea Ice polar stereographic north",
##         METHOD["Polar Stereographic (variant B)",
##             ID["EPSG",9829]],
##         PARAMETER["Latitude of standard parallel",70,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8832]],
##         PARAMETER["Longitude of origin",-45,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8833]],
##         PARAMETER["False easting",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",0,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["easting (X)",south,
##             MERIDIAN[45,
##                 ANGLEUNIT["degree",0.0174532925199433]],
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["northing (Y)",south,
##             MERIDIAN[135,
##                 ANGLEUNIT["degree",0.0174532925199433]],
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["Polar research."],
##         AREA["Northern hemisphere - north of 60°N onshore and offshore, including Arctic."],
##         BBOX[60,-180,90,180]],
##     ID["EPSG",3411]]
extent(seaice.jan1980)
## class      : Extent 
## xmin       : -3400000 
## xmax       : 3600000 
## ymin       : -4875000 
## ymax       : 5850000
#reproject to WGS84 to fit with basemaps
seaice.jan1980.WGS84 <- spTransform(seaice.jan1980, CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))

#check crs and extent
crs(seaice.jan1980.WGS84)
## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs 
## WKT2 2019 representation:
## GEOGCRS["unknown",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]],
##         ID["EPSG",6326]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433],
##         ID["EPSG",8901]],
##     CS[ellipsoidal,2],
##         AXIS["longitude",east,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]],
##         AXIS["latitude",north,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433,
##                 ID["EPSG",9122]]]]
extent(seaice.jan1980.WGS84)
## class      : Extent 
## xmin       : -180 
## xmax       : 179.7 
## ymin       : 38.01 
## ymax       : 83.67
#tidy data for working in ggmap and leaflet
broom::tidy(seaice.jan1980)
## Regions defined for each Polygons
## # A tibble: 4,470 × 7
##       long     lat order hole  piece group id   
##      <dbl>   <dbl> <int> <lgl> <fct> <fct> <chr>
##  1 1000000 5850000     1 FALSE 1     0.1   0    
##  2 1075000 5850000     2 FALSE 1     0.1   0    
##  3 1075000 5750000     3 FALSE 1     0.1   0    
##  4 1050000 5750000     4 FALSE 1     0.1   0    
##  5 1025000 5750000     5 FALSE 1     0.1   0    
##  6 1025000 5725000     6 FALSE 1     0.1   0    
##  7 1000000 5725000     7 FALSE 1     0.1   0    
##  8 1000000 5750000     8 FALSE 1     0.1   0    
##  9  975000 5750000     9 FALSE 1     0.1   0    
## 10  975000 5800000    10 FALSE 1     0.1   0    
## # … with 4,460 more rows

Make basemap

#manual bounding box
seaice.box <- matrix(c(-180, 180, 38.01, 83.67),
                     ncol  = 2,
                     byrow = TRUE)
colnames(seaice.box) <- c("min", "max")
rownames(seaice.box) <- c("x", "y")
seaice.box
##       min    max
## x -180.00 180.00
## y   38.01  83.67
#Get the base map (foundational layer)
seaice_base.map <- get_map(
                       location = seaice.box,
                       source   = "stamen",
                       maptype  = "watercolor",
                       crop = TRUE
                      )
## Source : http://tile.stamen.com/terrain/3/0/0.png
## Source : http://tile.stamen.com/terrain/3/1/0.png
## Source : http://tile.stamen.com/terrain/3/2/0.png
## Source : http://tile.stamen.com/terrain/3/3/0.png
## Source : http://tile.stamen.com/terrain/3/4/0.png
## Source : http://tile.stamen.com/terrain/3/5/0.png
## Source : http://tile.stamen.com/terrain/3/6/0.png
## Source : http://tile.stamen.com/terrain/3/7/0.png
## Source : http://tile.stamen.com/terrain/3/8/0.png
## Not Found (HTTP 404). Failed to aquire tile /terrain/3/8/0.png.
## Source : http://tile.stamen.com/terrain/3/0/1.png
## Source : http://tile.stamen.com/terrain/3/1/1.png
## Source : http://tile.stamen.com/terrain/3/2/1.png
## Source : http://tile.stamen.com/terrain/3/3/1.png
## Source : http://tile.stamen.com/terrain/3/4/1.png
## Source : http://tile.stamen.com/terrain/3/5/1.png
## Source : http://tile.stamen.com/terrain/3/6/1.png
## Source : http://tile.stamen.com/terrain/3/7/1.png
## Source : http://tile.stamen.com/terrain/3/8/1.png
## Not Found (HTTP 404). Failed to aquire tile /terrain/3/8/1.png.
## Source : http://tile.stamen.com/terrain/3/0/2.png
## Source : http://tile.stamen.com/terrain/3/1/2.png
## Source : http://tile.stamen.com/terrain/3/2/2.png
## Source : http://tile.stamen.com/terrain/3/3/2.png
## Source : http://tile.stamen.com/terrain/3/4/2.png
## Source : http://tile.stamen.com/terrain/3/5/2.png
## Source : http://tile.stamen.com/terrain/3/6/2.png
## Source : http://tile.stamen.com/terrain/3/7/2.png
## Source : http://tile.stamen.com/terrain/3/8/2.png
## Not Found (HTTP 404). Failed to aquire tile /terrain/3/8/2.png.
## Source : http://tile.stamen.com/terrain/3/0/3.png
## Source : http://tile.stamen.com/terrain/3/1/3.png
## Source : http://tile.stamen.com/terrain/3/2/3.png
## Source : http://tile.stamen.com/terrain/3/3/3.png
## Source : http://tile.stamen.com/terrain/3/4/3.png
## Source : http://tile.stamen.com/terrain/3/5/3.png
## Source : http://tile.stamen.com/terrain/3/6/3.png
## Source : http://tile.stamen.com/terrain/3/7/3.png
## Source : http://tile.stamen.com/terrain/3/8/3.png
## Not Found (HTTP 404). Failed to aquire tile /terrain/3/8/3.png.
ggmap(seaice_base.map) +
  geom_polygon(data = seaice.jan1980.WGS84, aes(x = long, y = lat)) +
  coord_map(projection = "orthographic")
## Regions defined for each Polygons
## Coordinate system already present. Adding new coordinate system, which will replace the existing one.

leaflet(data = seaice.jan1980.WGS84) %>% 
  addProviderTiles("CartoDB.Positron") %>% 
  addPolygons()